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Abstract 



The present work deals with the problem of extracting a sky map of a particular emission process that has been 
observed in an experiment with different detectors each of them having a different spectral response. This is the 
arena of the so-called "methods of component separation", especially in the field of microwave experiments, like 
the Cosmic Background Microwave (CMB) missions COBE, WMAP or, at present, the Planck surveyor mission. 
For the Internal Linear Combination methods, the difference with respect to other approaches is that it uses only the 
spectral behaviour of the sought component as an input. This idea has been applied to the case of CMB emission. 
Since this emission is itself the calibration emission in those maps, the problem is simplified. In this work, we derive 
the general expression for a generic spectral behaviour of the sought emission. We also apply the method to some 
of the common missions in the range of microwave and sub-mm emissions: Galactic dust, Synchotron emission, 
Sunyaev-Zel'dovich effect and as a check for consitency to CMB also. The data are simulations that resemble those 
performed presently by the Planck surveyor mission. Moreover, we will also show how it is possible to optimize the 
extraction of the chosen emission by minimizing the output noise and the bias in the extraction of the component. 
Therefore, we call this method MILCA: Maximum Internal Linear Component Analysis. 
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1. Introduction 

The Planck space mission 

(iThe PLANCK coUaborationh was launched on th 
14th Ma y 2009, is the thi rd space mission, aft er the 
COBE dSmoot et al. 1 192l: iBennett et al. 1996h and 
the WMAP (iHinshaw et al. 20031: iHinshaw et al. 20091 : 
iBennett et al. 2003h . The Planck mission is a typical 
CMB experiment, where the sky is observed at different 
frequencies. The Planck mission is observing the whole 
sky at 9 fequency channels, 3 of them at 30 GHz, 
44 GHz and 70 GHz constitute the so-called Low 
Erequency Instrument, LEI, and the other 6 channels, 
at 100 GHz, 143 GHz, 217 GHz, 353 GHz, 545 GHz 
and 857 GHz, constitute the so-called High Erequency 
Instrument, HEI. 

In the following, we will write the value of the emis- 
sion of a all sky map as: 

T = A.S-\-N (1) 

where T is vector containing all the observations at the 
different channels, A is a matrix, called "the mixing ma- 
trix". It has dimensions of ritXHs, where rit is the number 
of detectors, or elements in the vector T, and ^^ is the 
number of components, or elements in the vector S. Sis 
a vector whose components are the values of the com- 
ponent to be extracted at each of the detectors and A^ 
is the noise contribution. We use the notation A.B to 
denote the usual matrix product. 



In other component separation techniques, 
such as Maximum Likelihood techniques 



( Eriksen et al. 2008h or Maximum of Ent ropy method 
(IHobson et al. 19981: IStoliarov et al. 2002h one uses a 
priori knowledge about the component term and/or 
mixing matrix A. On the other hand, the ILC technique 
(iBennett et al. 20031) assumes a limited knowledge 
on the mixing matrix A: only one component have 
a known spectrum. In the usual case of the CMB 
emission, one readily has that the spectrum of 
the source is fiat, simply due to the fact that the 
data are calibrated with respect to the CMB itself. 
Einnaly techni ques like the Independent Component 
Analysis ICA (iMaino et al. 20021: ICardo so et al. 20081: 
iDelabrouille et al. 20031) . renders the matrix A as an 
unknown. 



We will show in this paper, step by step from ILC 
to MILCA, that it is possible to minimize the effect of 
the bias and noise by imposing additional conditions in 
the determination of the wanted component. The same 
idea is valid both for the case of a diffuse emission or a 
point- source one. 
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2. Generalisation of ILC for a non-flat 
spectrum emission: the noiseless case 

In this section, we assume that there is no noise term in 
the map-making or it can be disregarded at a first ap- 
proximation and let for the next one the case with the 
addition of noise. Then , one has 



T = AS, 



(2) 



where, we remind that T is a vector containing the ob- 
served values for each of the channels. A is the so-called 
mixing matrix and contains the frequency dependence 
of the diff'erent physical components. It is defined as: 



f Fi(y)Hj(y)dy 



(3) 



where Fi(y) is the frequency dependence of the i-th 
component and Hj(y) is the spectral response of the 
detector j. , T is the vector containing the observed 
values in the map for each of the channels. Finally, S is 
the vector containing each template map of the physical 
eff'ects and for the diff'erent observed frequencies 
(channels). 



2.1. Method 



Let us consider 

i 



(4) 
(5) 



where 5"^ is the component we want to extract. The 
vector w contains the set of weights that will be used to 
obtain that particular component by linear combination. 
Now we will describe how w is determined. 

In order to obtain w, we proceed as in the well- 
known ILC method for CMB emission but allow for a 
non-flat spectrum in the emission law of the component 
to be extracted. That is, w will be determined by miniz- 
ing the variance under the following constraint 



Z^- 



1, 



(6) 



where / is a vector containing the frequency depen- 
dence of the component to be extracted at each channel. 
The vector / is related to A by / = ^^.A, where ec is a 
column vector whose values are except for the "c"-th 
entry, associated to the component to extract, for which 
its value is 1 . 



The variance to be minimized is given by 
V(Sc) = w^.T.T^.w, 



(7) 



where V(S c) is the variance of the component extracted 
by a linear combinaition of observed channels. We solve 
Eq. d?]) by using Lagrange multipliers 



'VV(Sc)-A.Vg 



(8) 
(9) 



where A is some constant. 

The gradient of V(S c) and g are given by 

V.V(Sc) = 2.Ct.w, 
V.g = /, 



(10) 
(11) 



where Ct is the covariance matrix of the observed map. 
Eq. Q can be expressed in matricial notation 



2.Ct -f 
f 







(12) 



This system of equations can be solved by using a 
block inversion method, whose general expression can 
be written as 



^ ,_i _ / A-1 + A-i .B.{D - C.A-^ .B)-^ J 
^ - [ -(D-C.A'.B)KC.A' 



(13) 



C. A^ -A^ .B.(D - C.A-^ .B)-^ 
(D-C.A KB)^ 



With this in hand, the result for the w estimator is given 

by 

- cJ.f 

w = l^\ (14) 

fCj f 

The same results h as been obtain in parallel by 
(iRemazeilles et al. 201Qh . 

One can readily see that for a flat spectrum, as it 
is the case of the CMB emission, we obtain the same 
results as the well-k nown ILC method Eq[T5] given by 
(lEriksenet al. 20041) . i.e. 



w 



Cj .1 

rr 1 

l^.Cj .1 



(15) 



From Eq. ([T4l) we can derive the estimator of Sc 
which is given by 



^ f^'r 



fi^T-f 



.T 



(16) 



2.2. Properties of the noiseless estimator 

In this section, we develope the explicit expression of 
Eq. (O for the case of Eq. Q. For this case, we obtain 



Ct = A.Cs.A^ 



(17) 



where Cs is the covariance matrix associated to each 
component. 

Some precaution must be taken before inverting Ct- 
This matrix has dimension ntXnt but has rank min(ns, nt) 
by construction. 

There are three possible cases. 

- The case ns > nf. we have a lower number of detec- 
tors than the number of components. In this case, it 
is not guaranteed that a linear combination can pro- 
vide an isolation of the wanted component. Indeed, 
in this situation the ILC method is not efficient and 
can provide severe bias. 
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- The case Us = ^^ we have the same number of de- 
tectors and physical components. In this case A is a 
square matrix, so there is no ambiguity when invert- 
ing Ct- The procedure is straightforward. 

- The case Us < nf. in this case Ct has rank Us < Uf. 
Consequently, the matrix is singular and can not be 
inverted. However, one may take its pseudo-inverse 
as defined from well-known singular value decom- 
position (SVD) defined as Cj" = U.DKU^, where 
Z> is a diagonal matrix containing the singular val- 
ues of Ct and U is an orthogonal matrix obtained 
from singular value decomposition of Ct- With this 
in hanf, one can write Ct = U.D.U^. We remark 
that in this case D is uniquely defined, but U is not 
and has still rit - ris degrees of freedom showing the 
fact that there are indeed multiple linear combinai- 
son of T elements which extract the wanted compo- 
nent. 



In the sequel, we will only consider the cases ris < 



rif 



The results we obtain for the S c estimator is given by 

^ l^.CA.C^.A^)-! 



-.A.5 



, — 1 



f = A.ec 

a 

y-s )cc J*' 



(18) 

(19) 
(20) 
(21) 

(22) 

(23) 
(24) 



2.3. Application to simulated data 

In this section, we apply the expressions obtained so far 
to the case of dust emission. We do not apply it to other 
emission laws because we wish to show how the non- 
inclusion of noise bias the estimator and therefore the 
previous equations, Eqs. (l24l) . need to be modified to 
incorporate the noise term, and the case of dust emission 
has revealed to be the clearer case. 

In fig 121 we present the results of extracting a dust 
component from simulated PLANCK observations at 
100 GHz, 143 GHz and 217 GHz whith a simulated 
CMB (using synfas t with the best fit of the WMAP 
7 years results, see (iLarson et al. 20091) ). a dust model 
(FDS8, (1)) and a synchrotr on template (we use Haslam 
et al. map a t 408 MHz (iHaslam et al. 1982h in its 
Healpix form (ILAMBDA websiteh ). Each of these tem- 
plates are shown in Fig[TJ In order to simulate the emis- 
sion at 100 GHz, 143 GHz and 217 GHz, we assume a 
constant spectral index for the synchrotron of - 3.0 and 
a grey body at 17 K and with a constant spectral index 
of -h 2.0 for the case of dust. Furthermore, all the maps 
are finally degraded to the 100 GHz beam of 9.5 arcmin 
r solution. 



In fig 121 we show the importance of the bias due to 
spacial correlations that is present after applying the 
noiseless version of the ILC algorithm. Indeed, using 
a block inversion procedure for Cs it can be shown that 
the estimator of 5" ^ is unbiased only in the case where 
the component we want extract is uncorrected with 
the other components. This can be shown from Eqs[T4l 
There, one can see that if blocks B and C of the ma- 
trix are null, the corresponding blocks B and C of the 
inverted amtrix will be as well null. Therefore, (C^M 

will be null. If 5"^ is correlated with some component, 
a bias will result as shown at eq|24l This bias comes 
from the fact that the minimum of the variance is not 
an eigenvector of the component space. In any practical 
case Cs will be correlated foregrounds to some extend 
(but for the CMB), especially in the Galactic plane. Of 
course, a simple way to diminish the eff'ect of bias is in- 
troducing a mask that would remove most of the regions 
where correlation is higher, close to the Galactic plane. 
Yet, we will present in the following question that noise 
can be included in the algorithm improving the results 
on the bias. 

3. Noise induced bias 

Let us remember Eq. Q, 

T = A.S-^N. 

The vector A^ is a vector containing the noise component 
for each of the channels. In the noisy case Ct becomes 

Ct = A.Cs.A^ + Cn, (25) 

where C^ = N.N'^ is the noise covariance matrix. 

In this case, the matrix Ct is invertible. If we fur- 
ther assume as a first hypothesis that instrumetal noise 
among detectors is not correlated at first order C^ be- 
comes diagonal. Therefore, since C^ is a rit rank ma- 
trix, so it will be Ct- However, for Us < rit we can no 
longer describe Ct within a ris dimensional subspace. 
In this situation, the only way would be to increase the 
dimension of the components space to achieve a dimen- 
sion equal to the number of detectors. However, in this 
space Cs has no longer the physical sense previously 
defined and will unavoidably produce some extra bias 
on the S c estimator. Furthermore, Cn will contribute to 
the non diagonal terms in Cs resulting in a subsequent 
incrementation of the bias given by 



-1 



Sc = 



fCr 



T, 



rp 1 






(A.S -h A^). 



(26) 

(27) 



(28) 



In fig|3lwe show the bias due to the introduction of 
noise. The frequency bands are the same as in Fig 121 
In this case, the bias remains around 12% of the CMB, 
which is a huge bias. 

In any realistic situation, bias in the estimator can 
be produced by several diff'erent causes. For instance: 
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Figure 1. from left to right, simulated CMB whitout noise at a resolution of 9.5 arcmin, dust FDS8 model at 353 
GHz and Haslam et al. map for synchrotron at 408 MHz. 




Figure 2. Left : Obtained map with the ILC algorithm for dust whitout noise. It has a total RMS = 0.8. Right : 
Residual map between the ILC map and the dust imput map at 353 GHz showing the bias. The total RMS is now 
0.4 




Figure 3. Left: ILC obtained map for dust with noise. The total RMS is 1.02 
derived map. the total RMS is 0.07, and the RMS of noise is around 0.07 



Right: bias and noise in the ILC 



Statistical errors on the estimation of Ct due to the 
limited resolution of the experiment. 
Calibration uncertainties that will produce slightly 
different frequency spectra for the emissions differ- 
ent to the calibration one. This can be translated into 
an error associated to the mixing matrix A. Yet, in 
the case of high signal to noise ratio, ILC is very 
sensitive to calibration errors as the expected ones in 
the Planck mission, around 1% (iDick et al. 20091) . 
Difference between true A and the "a priori" spec- 
trum used for the separation. 
Instrumental noise as we have shown in this section. 



All those effects will contribute to bias on the Sc esti- 
mator. We try in the next section to minimize the bias 
on the estimator. 



4, Optimised estimator 

The variance associated to the Sc 
pixel is given by eql3Ql 



V(p) = < Sc(p).Sc(p) > 



V(p) = 



fCj 



■.T(p).T\p). 



T 

Ct f 



estimator for each 

(29) 
(30) 



The error in this technique can be minimized by using 
the condition number k which bounds the relative error. 
Its expression is 

'< = \h — II' (31) 

where Amax is the maximum eigenvalue of Ct and Amin 
is the minimal eigenvalue. In a case of low noise and 
where ris < rit, the Ct estimator will be close to a sin- 
gular matrix (not singular due to the statistic variance 
associated with the Ct estimation and with instrumental 
noise), which may result in a severe increase of noise. 
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4.1. Unbiased estimator 

All the issues considered in previous section intro- 
duce an additionnal bias on the estimator of Sc- 
A way to create an unbiased estimator is given by 
(IVio & Andreani 200 9). 



Ct - Ct - C\ 



(32) 



Eql32l provides an unbiased estimator in the case 
when the covariance matrix of the noise is known and 
which may be assumed in many practical cases. In this 
solution 5"^ is only biased by the correlation between 
Sc and the other components. The result is the same 
as obtained in Eql24l Yet, clearly, one does no longer 
minize the variance of noise. This is presented on FigH 
where we have now use now 2 more frequency bands 
(353 and 545 GHz). 

In Fig.|4]we show the result of applying the unbiased 
estimator. It turns out that when removing the noise, we 
are close to the case where Ct is nearly singular. And, 
as expected, there is a severe increase of the noise in the 
derived ILC map. In fact, there is still a 5% of the CMB 
emission remaining in those maps. 

4.2. Compromise between noise and bias 

In the rest of this work we will denote Ct by Ct- 

The eigenvalues of Ct which are not too small give 
us an idea of the number of Us that may be extracted 
with a reasonable unbiased and signal to noise. 

In the case ris < fit, the pseudo-inverse of Ct can 
be defined whith rit - ris degrees of freedom, as ex- 
plained in preceding sections. In practice we invert Ct 
in the subspace of extractible component. Hence, we 
will use those remaining degrees of freedom for min- 
imizing the noise RMS eqOll The results obtained for 
dust are shown in FiglH We remark, though, that this 
technique will introduce some additional bias. In fact 
any ILC and in particular the present MILCA algo- 
rithm will always have to find a compromise between 
the noise RMS minimisation and the bias removing. In 
our case, this translates into 



Vn = w .N.w 

g' = vy^ = 0, 



(33) 
(34) 



w:w here g' is a new constraint, and w has only rit - ris 
degrees of freedom. We remember that the minimiza- 
tion of V(Sc) is giving already ris constraints. In pratice, 
we have minimized g' using markov chains. We present 
the obtained MILCA map in fig IS) We obtain a quite 
well compromise between noise and residual bias. 

The Fig 15] shows that the residual bias is however 
important. Adding some more knowledge about com- 
ponent maps will improve the results as we detail in the 
next section. Here there is a 9% remaining CMB contri- 
bution. 



4.3. Improvement by acting on covariance matrix 
eigenvalues 

In the case of an unbiased estimator by removing noise 
covariance matrix, the eigenvalues of the covariance 
matrix correspond to the variance associated to the 
component (in the case where the components are 
statistically independent from each other). For usual 
ILC algorithm the component who dominates the 
signal in the map will have the best efficiency in the 
ILC procedure. Therefore, by artificially increasing 
an eigenvalue of one of the components, the MILCA 
algorithm will be minimize this component with a 
better efficiency. 

For instance, we compute the MILCA extraction in- 
creasing the eigenvalues associated to all components 
but dust to have a value 10 time greater than the one 
associated to dust. 

As shown in Fig. [6l we can also use this method 
to reduce the contribution of a component with an 
unknow spectrum. In the maps in Fig. [6l the remaining 
CMB is of 0.34%. This can be very helpful in the case 
the MILCA component extraction is to be applied in 
small patches. In this case, one can remove sources 
with diff'erent sort of spectra. This framework will be 
the matter of the remaining of this paper. We finish 
this section mentioning that by increasing some of the 
components, one necessarily degrades the capability 
of removing the other ones. Therefore, this procedure 
should not be carried in the case one has a biased 
estimator, because decreasing the noise would turn into 
an increase of the noise induced bias. 



5. Using additionnal (knowledge of the other 
components 

An additionnal way to remove the bias is using a "a 
priori" knowledge about A, as mentioned in preced- 
ing sections. In the case with ris < ^^ one has to use 
the pseudo inverse of Ct- We define / as a rectangu- 
lar matrix (rit^ric) with ric < ris, that will contain the 
frequency dependance of "a priori" known component. 
Incidentally, in this way, one also has the possibility of 
forcing an unwanted component to be null in the final 
linear combination. 

The new system of equations to solve is 



VV(Sc)-A^.Vg = 0, 



g=f.i 



ei. 



(35) 
(36) 
(37) 



with A a lie vector of Lagrange multipliers, 1 a vector 
of n, dimension containing only values 1 and ei a tic di- 
mension vector containing values 1 for the chosen com- 
ponent to be extracted, and for the rest. The new sys- 
tem of equations in matricial notation is 



2.Ct -f\lw 



(38) 
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Figure 4. Left: obtained map for dust unbiased estimator using the ILC algorithm. The total RMS is 9.56. Right: 
Residual map showing the bias. Total RMS is 9.54. The RMS of noise is 9.53 




Figure 5. Left: Obtained map with the MILCA algorithm for dust after minizing for the noise RMS. The total RMS 
is L02. Right: bias in MILCA map. The total RMS is 0.01 and the noise RMS = 4 lO'^ 




i^m'm^r^m^^^ 




Figure 6. Left: Obtained map using the MILCA algoritm for dust with noise. Total RMS is 1.02. Right: bias and 
noise from the MILCA map. The corresponding RMS are 5.10"^ and 3.10"^, respectively. 



Following the same procedure as in the previous 
section, we obtain an estimator for Sc given n Eq. [39l 
with the simple replacement of / by /. Consequently, 
the properties of this estimator are the same as the ones 
in the previous section, but have the remarkable advan- 
tage of having removed the bias due to the know un- 
wanted component. 



known component by this procedure eql40l 



((/^.' 



Ct ./)" 



'.r'c, 



V 



(39) 



Following a similar idea (adding extra-constrain on 
ILC like algorit hm), an other esti mator has been pro- 
posed by (Remazeilles et al. 20101) . 
Where now subscript c, refers to the "c"-th line of the 
matrix. The corresponding results are shown in FigUl 
One can see how there is no more a remaining residual 
from the CMB emission (assuming of course there are 
no calibration errors). Clearly, one can also extract other 



s = (f.c;\fr\r^.c;\T 



(40) 



6. The case of the patches approach. 
Application to the thermal SZ recovery 

In the preceding sections we have demonstrated that the 
MILCA method gives good results when appplied to a 
all sky map to recover diffuse emissions. The case of 
recovering sources can be indeed improved by using the 
MILCA algorithm in small patches around the source, 
which is a natural and easy to implement option, once 
the source candidates are kno wn from som e catalogue 
or method (as e.g. SExtractor (JBertin et al.l) ). 

In fact, this is the only left to proceed since the 
spectra of the sources covers many different values and 
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Figure 7. Left: Obtained map with the MILCA algorithm for dust whith a suppression of CMB emission. The total 
RMS is L02. Right: bias in MILCA map. The RMS is S.IO'^ and the noise RMS is 3A0-\ 



the number of components to extract cannot exceed the 
number of detectors. 

Instead of showing how the MILCA algorithm 
works for some sort of predefined sources with 
some known spectrum (as e.g., supernova remnants, 
Gigahertz peaked sources, dusty galaxies, ...) we will 
apply it to the more common interest case in CMB stud- 
ies of the thermal-SZ emission. 

To that end, we have performed a all-sky simula- 
tion, and have modified the thermal dust emission so 
that when selecting diff'erent patches we will also test 
the methods against variable dust index (between L4 
and 2.0 smooth at a resolution of 3 degrees) over the 
sky. 

In Fig. [8] we present the results of the extraction of 
the SZ signal for two diff'erent cases: if one applies the 
method to the whole sky and when applied on L7 de- 
grees square patches. We obtain lower noise results for 
the case of the L7 degrees square patches. The reason 
is again in the fact of the available number of indepen- 
dent observations (detectors) and the eff'ective number 
of components. In the case of real thermal dust emis- 
sion, the dust component cannot be described by only 
one spatially constant spectrum, which results in an im- 
possibility for any global linear algorithm to remove 
dust on the whole sky at once. However, as it is clear, 
when considering smaller regions, as the L7 degrees 
patches the dust spectral index will be in many cases 
practically constant and the algorithm will provide a 
better removal of this component. 

On the other hand, we remark that the small patches 
approach has obviously the drawback that the un- 
certainty of the covariance matrix estimation will be 
higher, and can yield to some artificial correlation be- 
tween an uncorrected component such as dust, CMB 
and thermal SZ. So, finally, one should find a compro- 
mise between the size of the patch and the limited num- 
ber of spectra that can be removed by an ILC like algo- 
rithm, as the MILCA algorithm is. Nevertheless, let us 
mention again, that we have provided the tools and the 
main ideas to be applied to each particular, real, situa- 
tion. 



7. Conclusion 

We have presented a method for extracting a chosen 
component from a map which is the result of the contri- 
butions of diff'erent components, as is the usual situation 
in the area of microwave and submm observations, like 
all CMB experiments. 

In a first step, we have used our knowledge on the 
frequency dependence on ric known component for con- 
straining the weights that allow to extract that compo- 
nent from the mixture by an, internal, linear combina- 
tion. 

In a second step, we minimize the variance of the 
chosen component obtaining ris - ric constraints on the 
weight of MILCA. 

Finally, we use the last rit - ris degrees of freedom 
for minimizing the noise covariance matrix. 
As a result, we have shown how we obtain a good com- 
promise between bias removal and noise RMS mini- 
mization. 

As for improving the present results, we would like 
to mention that the MILCA method can also be ap- 
plied on the Four ier space, in similar way as in e.g. 
(iTegmark et al. 20 03), In fact for improving results, the 
algorhitm can applied with some window both in real 
and Fourier space. Also, we would like to compare the 
results when using patches over all sky and when ap- 
plied directly to the all sky. 

For the extreme case where ric = ^s, our technique 
converges to the same solution that the one obtained by 
a maximum likelihood method (see Eq. (l4T1) for f = A 
and He = Us < fit). In this case, we can no longer con- 
strain the w matrix using the minimisation of the vari- 
ance of components. The last fit - fis degrees of freedom 
are constraint by minimizing the variance of noise, as it 
is done for the maximum of likelyhood method. 



S = {A^.C^ .Ay 



-^A^ 



.C^.T 



(41) 
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